test_mod1 <- lm(java_index_skd_bw_1 ~ fail_skd + AGE, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod2 <- lm(non_java_index_skd_bw_1 ~ fail_skd  + AGE, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod3 <- lm(regional_index_skd_bw_1 ~ fail_skd  + AGE, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod4 <- lm(religious_index_skd_bw_1 ~ fail_skd  + AGE, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod5 <- lm(corruption_index_skd_bw_1 ~ fail_skd  + AGE, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod6 <- lm(national_index_skd_bw_1 ~ fail_skd  + AGE, data = data_1 %>% filter(abs(forcing_skd) < 5))

observations <- c(nobs(test_mod1),nobs(test_mod2),nobs(test_mod3),nobs(test_mod4),nobs(test_mod5),nobs(test_mod6))

test_mod1 <- coeftest(test_mod1, vcov=vcovHC(test_mod1,type="HC0"))
test_mod2 <- coeftest(test_mod2, vcov=vcovHC(test_mod2,type="HC0"))
test_mod3 <- coeftest(test_mod3, vcov=vcovHC(test_mod3,type="HC0"))
test_mod4 <- coeftest(test_mod4, vcov=vcovHC(test_mod4,type="HC0"))
test_mod5 <- coeftest(test_mod5, vcov=vcovHC(test_mod5,type="HC0"))
test_mod6 <- coeftest(test_mod6, vcov=vcovHC(test_mod6,type="HC0"))

table <- list(test_mod1, test_mod2, test_mod3, test_mod4, test_mod5, test_mod6)

note_text <- paste(" Beta coefficients from OLS regression with age control. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measure are indexed values capturing (1) Javanese preferentialism among Javans and (2) among non-Javans, (3) regional preferentialism,
                   (4) religious resentment, (5) perceptions of corruption, (6) national identification.")

table = stargazer(table, 
                  type = 'latex', 
                  title = "The Effect of Basic Competence Examination (SKD) Failure, Age Control",
                  label = 'tab:testing_effect_age',
                  model.names = F,
                  model.numbers = T,
                  digits = 3,
                  column.separate = c(2, 1, 1, 1, 1),
                  column.labels = c("Java. Pref.", "Reg. Pref.", "Relg. Resent.", "Corrup. Percep.", "Natl. ID"),
                  
                  multicolumn = T,
                  dep.var.labels = NULL, 
                  add.lines = list(c("Subset", "Javan", "non-Javan", "---", "---", "---", "---"),
                                   c("Observations", observations),
                                   c('Bandwidth', rep(c('1\\%'), 6))),
                  covariate.labels = c("Failed Test", "Age"),
                  #star.cutoffs = c(0.05, 0.01),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')


write_latex(table[-c(10, 11, 12, 18, 21, 24, 31)], note_text, './_4_outputs/tables/table_a7.tex', .8)

